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Recent experimental work in the field of synthetic protocell biology has shown that prebiotic vesicles are able 
to 'steal' lipids from each other. This phenomenon is driven purely by asymmetries in the physical state or 
composition of the vesicle membranes, and, when lipid resource is limited, translates directly into 
competition amongst the vesicles. Such a scenario is interesting from an origins of life perspective because a 
rudimentary form of cell-level selection emerges. To sharpen intuition about possible mechanisms 
underlying this behaviour, experimental work must be complemented with theoretical modelling. The aim 
of this paper is to provide a coarse-grain mathematical model of protocell lipid competition. Our model is 
capable of reproducing, often quantitatively, results from core experimental papers that reported distinct 
types vesicle competition. Additionally, we make some predictions untested in the lab, and develop a general 
numerical method for quickly solving the equilibrium point of a model vesicle population. 

A fundamental problem in biology concerns the origins of an innovation that allowed the development of 
organisms in our biosphere, beyond complex chemical reaction networks: the emergence of cells 1,2 . Cells 
define a clear scale of organization and, given their spatially confined structure, they constitute efficient 
units where molecules can easily interact, coordinate their dynamical patterns and establish a new level of 
selection. Although it is often assumed that there was a transition from some type of 'less-organised' prebiotic 
chemistry (probably including catalytic cycles) to a cell-based living chemistry, little is yet known about the 
potential pathways that could be followed to cross it. Once in place, protocell assemblies would require available 
resources for their maintenance and, thus, would naturally get inserted in diverse competitive dynamics in which 
the main selective unit would be the whole protocellular system. In this context, aggregate-level evolution is the 
right scale of analysis to be considered. 

Different types of protocellular systems of diverse complexity have been studied from a theoretical stand- 
point 311 . In particular, by considering the coupling of a template carrying information with vesicle replication 
and metabolism, it has been shown that Darwinian selection is the expected outcome of competition in a 
protocellular world 12 . In a more simple scenario for autopoietic (i.e. self-producing) vesicles in a homeostatic 
regime, previous numerical simulations suggest that random fluctuations can also act as 'selection rules' for the 
more robust individuals 13 . Early pre-Darwinian stages in the development of biological organisms in which 
supramolecular systems could still be disconnected from information (i.e. closer to elementary forms of meta- 
bolism and strongly constrained by the molecular diversity of the available chemical repertoire) ought to be 
further explored. What type of competition and cooperation processes were at work in the chemical world leading 
to the emergence of early protocells? Processes able to favour asymmetries in the chemical composition of vesicles 
should be expected to play a relevant role in this context, defining the conditions under which protocellular 
assemblies could thrive. 

Recent laboratory experiments have actually demonstrated how differences in the composition or physical 
state of the vesicle membrane can drive competition for simple amphiphilic molecules (typically fatty acids), 
present in solution as free monomers at low concentration. First, Cheng and Luisi 14 observed competition 
between pure oleic acid and POPC vesicles, where each of these vesicle populations had different initial size 
distributions. In all the studied cases, the final size distribution was found near to the initial one of the POPC 
vesicles, suggesting that oleic acid molecules were rapidly absorbed by POPC aggregates. Then, Chen et al. 15 
reported competitive dynamics in a population of fatty acid vesicles, whereby vesicles that were osmotically 
swollen by an encapsulated cargo of RNA (or sucrose) stole lipids from their empty, osmotically relaxed counter- 
parts by virtue of absorbing monomers more quickly from the solution. They studied both oleic acid and POPC 
vesicles, but only in the former case was competition observed. This distinctive behaviour of fatty acid vesicles has 
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Table 1 | Vesicle competition model parameters. |a bl Area and volume of 1 20 nm diameter sphere. |c| 0.95Q 0 h reported in 24 . idl Close to 
80 fiM value reported for oleate vesicles 15 . |e| Bicine buffer concentration used experimentally 1517 . " See Supplementary Material for 
rationale. 16,1 Reported in 33 . lh) Considered similar to POPC head area reported in 33 . H Calculated from equation (1 3); same for spherical 



and deflated vesicles. » Calculated as: N = L°,+P U where L° =2S°, / {ga P + a L ), P fl =gL° l and g = p/(l - p) 
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Inside/ouside surface area of extruded model vesicles 1 " 1 
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Critical vesicle concentration for oleic acid 1 " 4 

Buffer concentration 1 "' 
OA lipid release constant 1 ' 1 
OA lipid uptake constant' 11 
OA lipid head area' 9 ' 
DOPA lipid head area" 1 ' 

OA monomer equilibrium concentration, for pure OA vesicle 1 ' 1 

OA monomer equilibrium concentration, for DOPA: OA (po = 0. 1 , d = 0) vesicle 1 ' 1 

Aggregation number for extruded pure OA vesicle 1 ' 1 

Aggregation number for extruded DOPA : OA (po = 0. 1 , a 1 = 0) vesicle 1 ' 1 
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been theoretically rationalised 16 by assuming that double chain phos- 
pholipids are taken up from solution by the vesicle membranes five 
orders of magnitude more slowly than single chain fatty acid 
molecules. 

More recent experimental work has turned attention to other 
possible selective advantages of protocells, such as phospholipid- 17 
and peptide- 18 driven competition amongst vesicles. Instead of mem- 
brane tension, the main factor for competition here is the different 
composition of the membrane; single-chain fatty acids are mixed 
with double chain amphiphiles or with a different type of surfactant 
molecule, like sufficiently hydrophobic peptides. In the case of phos- 
pholipid-driven competition, oleic acid vesicles endowed with a 
membrane fraction of phospholipid are observed to take fatty acid 
molecules from phospholipid-deficient neighbours, who shrink, 
whilst the former grow and keep their potential for reproduction. 

In this paper we develop a mathematical model of a competing 
population of vesicles, with the motivation to explore and test pos- 
sible mechanisms underlying lipid competition phenomena. The 
model is based at the coarse-grain level of lipid kinetics, following 
the approach of Mavelli and Ruiz-Mirazo 16 . Using physically realistic 
parameters such as lipid molecule sizes, vesicle aggregation numbers 
and critical vesicle concentrations (CVC) as detailed in Table 1, we 
are able to qualitatively and often quantitatively reproduce results 
from two key experimental papers describing phospholipid-driven 17 
and osmotically-driven 15 competition. 

In the model, a vesicle in a population absorbs and releases amphi- 
philes to and from its membrane at rates that depend on the current 
physical properties of that particular vesicle (such as membrane 
composition or extent of osmotic tension). To take account of phos- 
pholipid-driven competition, we build into the lipid kinetics two 
basic physical mechanisms, which have been postulated in the lit- 
erature to underlie asymmetric growth dynamics in this context: the 
indirect effect and the direct effect, as will be named in this work. The 
first one refers to the decrease of amphiphile release processes simply 
due to the fact that other surfactant molecules are present in the 
membrane, and the second to the immediate influence that these 
surfactant molecules could have on the amphiphiles (see Fig. 1). 

More generally, this work forms part of our endeavour to try to 
develop a formalism that grasps the lipid kinetics involved in vesicle 
self-assembly under controlled conditions (pH, temperature, etc.). In 
contrast to the kinetics of chemical reaction networks which have 
been extensively modelled by the Mass Action Kinetics (MAK) and 
Stochastic frameworks 19,20 , membrane lipid kinetics have been lar- 
gely under- explored in the literature, due to the inherent complexity 



of supramolecular structures. Nevertheless, models coupling 
together membrane and metabolism kinetics will be a crucial corner- 
stone in order to build a systems understanding of the dynamic 
properties and organization of protocells, ultimately biological cells 
as well. 




Figure 1 | Two mechanisms of phospholipid-driven growth, (a) Indirect 
effect, whereby the presence of phospholipid in a vesicle membrane drives 
growth simply through a geometric asymmetry: only the lipid section of 
the bilayer (grey) is able to release lipid (orange arrows) whereas the whole 
of the bilayer surface (made of lipids and phospholipids) is able to absorb 
lipid monomer (green arrows). Phospholipid fraction is pictured as one 
continuous block to highlight the principle only. The indirect effect can be 
created also by non-lipid surfactant molecules (e.g. peptides) residing long 
enough in the membrane to increase surface absorption area, (b) Direct 
effect, whereby the acyl tails of the phospholipids have high affinity for 
packing closer to each other and increasing bilayer order, thus making the 
exit of the simple lipids more difficult. The direct effect is specific to the 
molecular structure of phospholipids. 
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The paper is organised as follows. The remainder of the introduc- 
tion serves to both introduce our kinetic model in detail, and to 
perform a mean-field analysis of it. This analysis gives insight into 
why we should expect phospholipid-driven competition to result in a 
simplified version of our model. Then, the Results section sum- 
marises how well numeric simulations of the full kinetic model are 
able to reproduce experimental results and observations, including 
also some predictions for still untested protocell competition scen- 
arios. In the Discussion section, we comment on some assumptions 
and other aspects of our approach and conclude the study. The 
Methods section at the end of the paper describes a fast numerical 
method for solving the final equilibrium state of the full model. This 
method was essential in producing the results figures in the paper. 
The Supplementary Material (online) explains aspects in more detail, 
including justification for some modelling choices and the vesicle 
mixing procedure assumed in order to compare our model with 
experimental observations. 

Theoretical model of vesicle competition. The competition model 
(Fig. 2) involves a set of n vesicles V = {Vi, . . . ,V„}, each one 

characterized by a quintuple Vj = \Qj,L^,P^,L{,BlA and all 

embedded in a finite volume environment £ defined by a triple 
(fl e , L e , B e ). 

Each competing vesicle Vj consists of a unilamellar (single bilayer) 
membrane of up to two different lipid types: single-chain fatty acid 
lipids V (e.g. oleic acid, OA), and possibly a fixed number of double- 
chain phospholipids P (e.g. di-oleoyl-phosphatidic acid, DOPA). 
Membrane thickness is considered negligible, and surface area of a 

vesicle, referred to as S^= — (l^Xl + P^-p^ , is the water-exposed 

area of either of the bilayer leaflets. The L type lipids in the bilayer 
continuously exchange with the vesicle internal water pool and £, 
whereas the phospholipids P are considered approximately station- 
ary in the bilayer due to their comparatively slow exchange rate, in 
agreement with previous reported work using POPC vesicles 16 . The 
internal water pool of each vesicle is considered a well-mixed chem- 
ical domain of volume Slj and hosts V c lipid monomers and also Bj. 
buffer species. Buffer species cannot permeate the bilayer but provide 
osmotic stability and they are also present in £ with constant number 
B, 

Vesicles compete with each other by consequence of uptaking/ 
releasing fatty acid monomers L from/to £ , which is a common 
limited resource. The initial system of vesicles is taken to be the result 
of mixing different vesicle populations, and is a closed system in a 
non-equilibrium state. The system equilibrates to a final state follow- 
ing the dynamics described below, with some vesicles growing bigger 
in surface at the expense of others, which shrink. We ignore spatial 
correlations and the possibility of direct vesicle-vesicle interactions, 
and assume a well-mixed set of vesicles. 

More precisely, each vesicle Vj is considered to release lipids to 
both aqueous phases (at each side of the bilayer) at the equal rate of 

= k 0U tlJ ll r\j)A, and absorb lipids from each phase at rate 

A'" = ki„8 ll [L]u(<f>]), where [L] is the molar concentration of lipid 
monomer in the respective phase. Functions r and u are defined later. 

The uptake and release kinetics are symmetric on each side of the 
bilayer, which means that the lipid monomer concentration inside 
and outside each vesicle will be equal [Lf c = [L] e = [L] at equilibrium. 
Flip-flop of the fatty acid L between membrane leaflets is considered 
very fast with respect to its uptake and release rates, and thus a bilayer 
is modelled as a single oily phase; this simplification is supported by 
experimental work from Hamilton's lab 21,22 . Conceivably, leaflet 
asymmetries could be created by the fact that the flip-flop of depro- 
tonated and protonated fatty acid molecules is not the same 23 . 



However, such effects are considered of secondary importance and 
are disregarded in the present work. 

Explaining the choice of L release kinetics, each fatty acid in a pure 
L membrane is considered to have a uniform probability per unit 
time k out of disassociating from the membrane 16 , while function r has 
been introduced in this work to take into account the direct effect. 
This function (0 £ r(p) £ 1) modifies the fatty acid release prob- 
ability, based on the current molecular fraction of phospholipid P in a 
p 

membrane p = — . It is monotonically decreasing with increas- 

P)i "1" hp 

ing p, meaning that increasing phospholipid fraction generally 
decreases bilayer fluidity, slowing down the rate of L release from 
the membrane 17 . In a first approximation, r was assumed linear: 

r[p) = l-dp (1) 

where parameter 0 £ d £ 1 tunes how the lipid release rate is affected 
by phospholipid content ( 1 being maximally affected and 0 being not 
ataU). 

Conversely, lipid uptake kinetics reflect that the probability of 
uptaking a lipid L to the membrane is proportional to the density 
of lipid monomer in the immediate vicinity of the respective bilayer 
surface (i.e. the concentration of lipid in the surrounding medium), 
the area of surface available for absorption S A , and function u, based 

on the dimensionless reduced surface <t> = S,, I \/367tn 2 . The reduced 

surface encodes the surface area to volume ratio of a vesicle, when the 
latter internal volume is considered as a sphere: Q> = 1 denotes a 
vesicle perfectly spherical in shape, whereas €>< 1 or <P > 1 indicates 
a vesicle in osmotic tension or deflated respectively. Taking this into 
consideration, we define u as the following conditional function 

u$= VK<t> ' 2 
[1, <D>1 

to denote that lipid uptake is only increased when the the bilayer is 
stressed 24 . Flaccid vesicles do not have extra enhancement of lipid 
uptake rate. Rationale for this function originated in the theoretical 
modelling of osmotically-driven competition dynamics between 
fatty acid vesicles 15,16 , and additional justification is provided in the 
Supplementary Material. 

The indirect effect is manifest as a systems property of the model, 
rather than in any particular function. When a vesicle membrane 
contains phospholipids (or other surfactant species like hydrophobic 
peptides), the P fl molecules add a contribution to the surface, increas- 
ing the L uptake rate, whereas the L release rate remains unaffected 
by their presence. 

Uptake and release kinetic constants k in and k out are set by two 
criteria. The first criterion is that pure fatty acid model vesicles (made 
solely of L), either spherical or deflated, must be in equilibrium when 
the fatty acid monomer concentration inside and outside the vesicle 
is the CVC for that amphiphilic compound (e.g. oleic acid). The 
second criterion is that the model dynamics must reproduce, with 
lowest RMS error, the experimental time courses reported by Chen et 
al. 15 for surface changes in osmotic competition. The second criterion 
narrows the possible {k in , k out } pairs (see Supplementary Material). 
For mixed membrane vesicles containing both L and P lipids, we 
assume that the lipid kinetics equations define what lipid monomer 
concentration inside and outside the vesicle [L] eq is necessary to keep 
the mixed membrane vesicle in equilibrium (however, in reality, the 
CVC of mixed lipid solutions is not a trivial matter 25 ). 

For the purpose of lipid competition, £ has a fixed volume of Q e 
litres. Each vesicle Vj has, in principle, a variable internal water 
volume of Qj = Q e (IJ c +S l J /(L e + B e ) litres. This volume value is 
based on the assumption that water permeates the membrane extre- 
mely rapidly, and ensures that the interior of each vesicle is isotonic 
with respect to £ at all times. However, since in real fatty acid vesicle 
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Figure 2 | Kinetic model of vesicle competition. ( a) Our model approach considers as a starting point a population of vesicles (of generally heterogeneous 
sizes and membrane compositions) in a well-mixed environment, (b) Each vesicle has a membrane composed of simple single chain lipids L, e.g. oleic acid 
(OA), and (c) sometimes more complex double chain phospholipids P, e.g. dioleoyl-phosphatidic acid (DOPA). (d) outlines the kinetic interactions 
between vesicles. Here two vesicles are displayed (bilayer cross sections not to scale). Vesicle 1, on the left hand side, has a mixed membrane with 
approximately 10 mol% phospholipids P (black) and the remainder single chain lipids L (grey). Vesicle 2 consists purely of simple lipids L. In the ensuing 
competition, phospholipid-laden vesicle 1 will grow at the expense of vesicle 2, which will shrink. 



solutions the concentration of extra compounds like counter ions, 
pH buffer, etc., is much higher than the concentration of free fatty 
acid monomers, we can reasonably assume that the aqueous volume 
of each vesicle is approximately constant at f2j=f2 e (Bj./£ e ). Thus, 
vesicle volume is largely determined by the number of buffer mole- 
cules a vesicle has trapped inside its internal water pool, with L flux to 
and from the water pool having marginal osmotic effects. 

To summarise, the state of the vesicle system is captured by enu- 
merating the number of lipids in each of the aqueous pools inside the 
vesicles, and in each of the vesicle membranes. The ODE system that 
describes the time behaviour of the entire vesicle population consists 
of 2n equations, two for each vesicle: 

^ = k 0Ut L{T - k in ^ [L]{.u((D ; ) (3) 
dV / \ 

-£ = -2k 0Ut L\,r[ P) ) +k*S>M+WM*}) W 

At the same time, the total number of lipids in the system L t is a 
conserved quantity set by the initial condition of mixing (see 
Supplementary Material), always equal to the number of lipid mono- 
mers in the environment L e , plus the number of lipids composing the 
vesicles: 

^+^(4+4)-^=° ( 5 ) 

Therefore, L e can be deduced from equation (5) once all L c and L ;I 
have been calculated at time t. Values of model parameters are given 
in Table 1. 



Mean field approximation. In the first instance, before performing 
any numeric simulations, why should we expect phospholipid 
fraction and surface growth to be correlated in the vesicle 
competition model? To answer this question, we can make a mean 
field approximation. This approach considers a reduced scenario 
where many details associated to the full model are ignored in 
order to keep only the logic of the problem (Fig. 3). 

The first simplification will be to ignore the internal structure of 
the vesicles, describing them instead as coarse-grained 'aggregates', 
denoted by pairs Vj = {Lj,Pj) > which contain just lipids and phos- 
pholipids. This step can be considered justified on the grounds that, 
at equilibrium, the amount of lipid monomer residing in the vesicle 
water pools (which typically have tiny volumes, around 1 quintil- 
lionth of a litre) is marginal as compared to the lipid composing the 
vesicle membranes. Since the internal structure or topology of the 
vesicles is disregarded, it actually amounts to treating them as elon- 
gated micelles or flat bilayers. 

The second simplification involves reducing the lipid uptake and 
release equations to their most basic form, independent of mem- 
brane tension (u(€>) = 1) and independent of membrane phospho- 
lipid fraction (r(p) = 1) respectively. Thus, the ODE system reduces 
to n simplified equations, where for each aggregate: 

= -k 0Ut Lj + ^k in (LjO. L + Pja. P )[L] e (6) 

Under these conditions, at equilibrium, the molar lipid concentra- 
tion in the environment [L] e = [L] eq is related to the number of lipids 
and phospholipids in an aggregate by the following function: 
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Figure 3 | Meanfield model of vesicle population dynamics. Some 
analytical treatment of the model is possible if vesicles lose their internal 
structure, and are just considered to exchange fatty acid with the external 
solution following simplified kinetic rate equations. 

For a fixed number of phospholipids Pj > 0, the mapping f:Lt—> 
[L] eq can be verified to be one-to-one, meaning that each aggregate is 
in equilibrium at only one specific outside lipid concentration, 
dependent on the number of lipids Lj it contains. Thus, no multiple 
equilibria of the population are allowed from this type of aggregate 
dynamics. 

Now consider two arbitrarily chosen aggregates i and j in the 
population of n aggregates, which are competing for lipid. Their 
ODEs, when written as: 



dL, 
lit 



dLj 
~dt 



= -koutLi + iiiLiOCL+PiUpjlLt- 



t Lj + ii(Lja. L +Pja P ) I ( - ^ L 



where n = k jn /2N A Q. e , are reminiscent of the Lotka-Volterra com- 
petition equations associated to species sharing and competing for a 
common set of resources 26 . If we look for the equilibrium solutions of 
the previous system, using dLJdt = dLj/dt = 0, we obtain 



Litx L + PiiXp Li 
Lj0i L + Pj0ip L 



(8) 



which leads to the following proportionality relation at equilibrium: 



(9) 



meaning that the final equilibrium vesicle sizes will be correlated with 
their respective numbers of phospholipids. Unless P, = P ; - one of the 
vesicles will be larger and the second smaller. For each pair (P ; , Pj) 
with Pj Pj a single solution is found. 

When functions u and/or r are not constant, unless they have a 
trivial form, it is generally not possible to show analytically what 
shape the correlation between phospholipid fraction and surface 
growth will take. However, in the Methods section at the end of 
the paper, we develop a fast numerical way to find the equilibrium 
configuration of the fully-fledged vesicle population model, with 
vesicles recovering their internal structure. As compared to numer- 
ically integrating the ODE set, the method provides the extra advan- 
tages of (i) being faster and thus scaling better for large vesicle 
populations and (ii) being able to calculate competition 'tipping 
points' (i.e. critical points that mark the transition between growing 
and shrinking) directly. 

In the following Results section, our fast procedure was used to 
perform accurate vesicle stoichiometry calculations, whilst simula- 
tions of the model dynamics were carried out using small popula- 
tions of vesicles and deterministically integrating the ODE set. 



Results 

Two competing populations: comparison with experimental 
results. Figure 4 compares predictions made by our kinetic model 
against experimentally reported surface growth of vesicles (assessed 
by a Forster resonance energy transfer assay, FRET) in phospholipid- 
driven 17 and osmotically-driven 15 competition. Two scenarios are 
forecast by our model: one where vesicles are spherical when 
extruded, and one where they are deflated by 5% when extruded, 
as generally observed experimentally 24,27 . The Supplementary 
Material details the vesicle population mixing procedure used to 
initialise our theoretical model, adjusting it to realistic experimen- 
tal conditions. 

Top figures 4a and 4b show the dynamics of surface area change in 
phospholipid-driven competition. Figure 4a details, in real time, the 
relative surface area of a tracked (surface area followed by fluor- 
escence probe) population of DOPA:OA (p 0 = 0.1) vesicles, when 
this population is mixed 1 : 1 with either pure OA vesicles, similar 
DOPA:OA (p 0 = 0.1) vesicles or simply buffer. In Fig. 4b, the 
tracked population is instead pure OA vesicles, which are mixed 
1 : 1 with the same three options outlined above. 

Whether starting with initially spherical vesicles, or vesicles 
deflated by 5%, execution of our lipid kinetics model correctly pre- 
dicts that when mixed 1:1, DOPA : OA vesicles steal lipid and grow 
(rising lines, Fig. 4a) at the expense of the pure OA vesicles, which 
shrink (falling lines, Fig. 4b). In this case, there is also fairly good 
quantitative agreement with the experimentally observed time 
courses (RMS error given as Supplementary Table S2). For the other 
cases, the kinetic model correctly predicts approximately no surface 
area change (no competition) when similar populations are mixed, 
or when a population is mixed with buffer. 

Middle figures 4c and 4d show phospholipid-driven competition 
from a different angle: that of vesicle stoichiometry. Stoichiometry 
explores the final equilibrium size of vesicles in a tracked population, 
when this population is mixed with a different population containing 
approximately R times as many vesicles. In this approach, the trend 
of final equilibrium surface area size versus mixing ratio is explored, 
rather than the dynamics on the way to equilibrium. Figure 4c details 
final surface area of a tracked population of DOPA : OA (p 0 = 0.1) 
vesicles, when this population is mixed 1 : R with a population of pure 
OA vesicles. Figure 4d details the opposite scenario, whereby the 
tracked population is OA vesicles, mixed 1 : R with DOPA : OA vesi- 
cles. The R = 1 cases in Figs 4c and 4d correspond to the surface sizes 
reached in the limit of time in Figs 4a and 4b, respectively. 

Calculating competition equilibrium by means of the fast com- 
putation approach outlined in the Methods section, we were able to 
verify that our model exhibits continual growth of DOPA : OA (p 0 = 
0.1) vesicles as more OA vesicles are added (Fig. 4c). If vesicles started 
at 5% deflation, the model matched the experimental data points 
even more closely. In the opposite scenario, we also verified that 
the model shows the same distinctive plateau in the shrinkage of 
pure OA vesicles as more DOPA : OA (p 0 = 0.1) vesicles are added 
(Fig. 4d). In the case of the latter figure, notably the indirect effect 
alone is sufficient to reproduce experimental results. 

Importantly, the general outcome of phospholipid-driven com- 
petition in our model is for DOPA : OA mixed vesicles to uptake 
lipid, grow in surface and to finish at high CD > 1 values (excess 
surface, flaccid), whereas pure OA vesicles lose lipid, suffer reduced 
surface, and finish at (D< 1 values (osmotically tense, spherical). This 
is observed experimentally, and indeed provides the basis for the 
conjecture that phospholipid-laden vesicles are more likely to divide 
spontaneously when gentle external shearing forces are applied 17 . 

Moving to osmotically-driven competition, Fig. 4e shows simu- 
lation of a swelled population of vesicles competing with an initially 
isotonic (non-swelled) population. Simulation outcomes match 
quite well the experimental best-fit time courses, in particular for 
the growth of the swelled vesicles (less accurately for the shrinkage of 
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Figure 4 | Comparison between kinetic model predictions and experimental results. In all plots, model vesicles extruded with 120 nm diameter surface 
area, either being spherical (€> = 1, green lines) or deflated by 5% (S> = 1.0348, red lines). Solid lines denote indirect effect only (d = 0) and dotted lines 
denote maximal direct effect (d = 1 ) present in DOPA : OA vesicles. Top plots show dynamics of phospholipid-driven competition. Experimental data 
points from Budin & Szostak 17 (Figs. 1A, IB therein) reproduced in background as coloured dots, (a) Surface change of model DOPA : OA (p 0 = 0.1) 
vesicles over time when mixed 1 : 1 with pure OA vesicles, with similar DOPA : OA (p 0 = 0. 1 ) vesicles, and with buffer. Model outcomes when mixing with 
similar DOPA : OA vesicles shown as horizontal grey line, and when mixing with buffer, as blue lines, (b) Surface change of pure OA vesicles over time 
when mixed 1 : 1 with the same three options as in (a). Model outcome when mixing with OA vesicles or with buffer shown as horizontal blue line. Middle 
plots show vesicle stoichiometry effects in phospholipid-driven competition. Supplementary Material defines our interpretation of vesicle mix ratio R in 
detail, (c) Continued average surface growth in fixed population of model DOPA : OA vesicles as more OA vesicles added at increasing mix ratio Rand (d) 
plateau of average surface shrinkage in fixed population of OA vesicles as more DOPA : OA vesicles added at increasing mix ratio R. Black markers with 
error bars reproduce experimental data points from Budin & Szostak 17 (Figs. 1C, ID therein). Bottom plots show osmotically-driven competition results, 
(e) Growth dynamics of model swelled OA vesicles and shrinkage of isotonic OA vesicles compared against best-fit exponential decay curves (dotted blue 
lines) to experimental data points from Chen et al. 15 (Figs. ID, IB therein), (f) Stoichiometry effects in osmotically-driven competition. Shrinkage of OA 
vesicle surface reaches a plateau as more swelled vesicles are added at increasing mix ratio R (note log scale). Black line and markers reproduce 
experimental results from Chen et al. 15 (Fig. 2A therein). Minimum 3) reached by model OA vesicles in (a)-(d) is 0.7692, in (e) is 0.7046. 
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the non-swelled vesicles). In any case, it must be noted that the 
original experimental data (yellow data points) has considerable 
variance. Then, Fig. 4f shows that the kinetics model qualitatively 
reproduces the stoichiometric observation whereby adding more 
swelled vesicles to a population of initially non-swelled vesicles will 
cause the shrinkage of the non-swelled vesicles to plateau, rather than 
to continue (note the logarithmic scale of Fig. 4f). Again, model 
outcomes are improved if vesicles start at 5% deflation. 

The general outcome of osmotically-driven competition in our 
model is for all vesicles to finish with different surface sizes (as for 
phospholipid-driven competition), but now, all vesicles also share 
the same €> < 1 value, indicating equal osmotic stress. This residual 
osmotic stress is also observed experimentally and stands as the main 
criticism of the osmotically-driven competition scenario. In order to 
divide, swelled vesicles would have to overcome a stronger energetic 
barrier, changing their stressed membrane state into one ready for 
fission, making this an improbable route to spontaneous vesicle 
reproduction 18 . 

Our kinetic model can also be used to make predictions or to find 
competition 'tipping points' in the more general scenario where 
completely heterogeneous populations of phospholipid-laden and/ 
or osmotically swollen vesicles compete for lipid (Figs. 5 and 6), even 
if some of these experiments have not been realised in the lab yet. 

Competition tipping points in diverse populations. Figure 5a shows 
that within a hypothetical population of model phospholipid-laden 
vesicles, where each vesicle has a randomly assigned phospholipid 
fraction in the membrane between 0 and 100%, the critical DOPA 
fraction needed for growth (tipping point), in this case, is just over 
58%. 

Figure 5b compares different heterogeneous populations compet- 
ing for phospholipid, and reveals an important observation: competi- 
tion is always context dependent. That is to say, a certain amount of 
membrane phospholipid does not guarantee a certain final surface 
area. Rather, final surface depends on the boundary conditions of the 
competition event (that is, the parameters influencing the solution of 
equation (15) in the later Methods section), which includes the num- 
ber and composition of competitor vesicles present. For example, 
population (i) in Fig. 5b has vesicles with low DOPA fraction as 
compared to vesicles in population (iv), yet in some cases the vesicles 
in the former population have larger final surface growth than vesi- 
cles in the latter. This concurs with the experimental observation that 
even small differences in phospholipid content should drive 
growth 17 . 

The dotted black lines in Figs 5a and 5b are the same competition 
events run when the direct effect is present, and maximally enabled 
(d = 1). The extent to which the direct effect affects vesicle growth 
must be made on a case by case basis, as it depends on the specifics of 
the competition event. For example, the direct effect has marginal 
influence on vesicle growth trends in the population shown in Fig. 5b 
(iii), but is more relevant in population (ii). The Supplementary 
material contains a recalculation of both Fig. 5a and 5b, if we further 
take into account the realistic constraint that vesicles will burst when 
osmotic tension exceeds a critical limit (<P < 0.7). 

Figure 5c shows that in a heterogeneous population where pure 
OA model vesicles are swelled to differing extents, vesicles with low 
initial (D values take lipid from those with higher (less swelled) <£ 
values, with the tipping point between growing and shrinking at 
(DjJ"' =0.85. As a last remark, orange crosses marked on Figs 5a 
and 5c show that full deterministic simulations of the model (run 
all the way to equilibrium) agree with and thus validate the 'Fast 
Computation of Competition Equilibrium' procedure outlined in 
the Methods section. 

Theoretical predictions beyond current experimental results. 

Finally, we were able to explore more widely some of the 



parameter space for phospholipid-driven and osmotically-driven 
competition, using our model to make some predictions. Figure 6a 
shows the stoichiometry results of phospholipid-driven competition 
in this wider context. A population of DOPA : OA (p 0 = 0.1) vesicles 
is mixed with a second population, but the phospholipid content of 
the second population, as well as the mixing ratio R, are varied. 
Taking a slice through the surface labelled 'popl' when p^" 1 ' 2 — 0 
shows the result reported as the solid red line in Fig. 4c. Figure 6b 
explores the stoichiometry of osmotically-driven competition in a 
similar way to phospholipid-driven competition. A fixed population 
of swelled vesicles is mixed with a second population, where the 
degree of swelling in the second population, as well as the mixing 
ratio R are varied. To conclude these predictions, Fig. 6c shows the 
effects of osmotically-driven versus phospholipid-driven competi- 
tion, still a completely unreported scenario in the experimental 
literature, whereby a population of swelled pure oleate vesicles 
competes for lipid with a population of DOPA : OA vesicles. The 
swelled oleate vesicles are able to steal lipid from the DOPA : OA 
vesicles, when the former have a high degree of swelling and the latter 
have a low DOPA fraction; otherwise, the DOPA : OA vesicles 
prosper in the competition. 

Vesicle bursting is an important consideration in Fig. 6. Com- 
petition predictions in Fig. 6a are only strictly valid when the popu- 
lation 1 surface is above the red box lines. Below these lines, vesicles 
in population 1 have excessive osmotic pressure (<1> < 0.7) and would 
likely burst, altering competition outcomes for population 2. 
Likewise, the population 2 surface in Fig. 6c is only drawn for values 
where oleate vesicles in that population have final (D > 0.7. Outside 
the extent of the population 2 surface, competition outcomes for 
population 1 should be treated with caution, as not all population 
2 vesicles will be intact. 

Discussion 

In this work we have presented a theoretical model of the transfer 
kinetics of single chain fatty acids between competing vesicles. We 
have shown that data coming from controlled laboratory experi- 
ments on phospholipid-driven competition and osmotically-driven 
competition can be reproduced fairly well by a set of physically-based 
rate equations describing the uptake and release of fatty acids for 
each vesicle. Furthermore, we have been able to predict the outcome 
of several yet-to-be-performed experiments. Thus, it is time to recap- 
itulate, considering possible limitations of our approach, clarifying 
several points that remain open, and giving a more general perspec- 
tive on the problem addressed. 

The main assumption we made when modelling phospholipid- 
driven competition is that the phospholipids are not released by 
vesicle membranes at the timescale of fatty acid transfer between 
the supramolecular structure (i.e., the closed membrane bilayer) 
and the aqueous solution (both inwards and outwards). This 
assumption is well founded on experimental evidence 14,15 and a pre- 
vious theoretical model 16 . Likewise, the assumption we make with 
osmotically-driven competition is that the extra buffer present inside 
the vesicles, swelling them, permeates very slowly through the bilayer 
membrane. In reality, the off-rate of a lipid molecule from a bilayer is 
inversely correlated with the number of carbon atoms in the acyl 
chain of the lipid concerned, and phospholipids do have a small non- 
zero transfer rate (with a half time from hours to days 16,17 ). If the 
much slower phospholipid transfer was included in our model, the 
equilibrium reached in the limit of time would always be that of a 
completely homogeneous population. This is because the P phos- 
pholipid would redistribute amongst the vesicles until all were equi- 
librated with the same phospholipid monomer concentration in 
solution [P] cq , which is trivially when all vesicles have the same 
fraction of membrane phospholipids P„. With no remaining asym- 
metries in P„ fraction to drive competition, all vesicles would finish 
with the same lipid composition and same surface size. The initial 
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Figure 5 | Lipid competition tipping points, (a) Phospholipid competition between 30 model phospholipid-laden vesicles each with a different DOPA 
fraction randomly assigned over the uniform interval 0 < p 0 < 1 and initially spherical, 120 nm diameter. Depending on initial DOPA fraction, each 
vesicle starts at a point on the horizontal blue line, and grows (green arrows) or shrinks (red arrows) to a point on the black line. The form of the black line 
is specific to this particular competing population, and is computed by equation (15). Competition 'tipping point' is shown by blue circle: any vesicle 
with p 0 > Pq"' = 0.584 gains lipids from its competitors. The solid black line shows relative growths when only the indirect effect exists (d = 0); for 
comparison, the dashed black line shows relative growths when the direct effect also maximally present ( d = 1 ) . Initial fatty acid concentration inside and 
outside model vesicles was [I] = 5.0 X 10~ 5 M. (b) Phospholipid competition in four unique populations of 30 model vesicles, with DOPA fraction 
randomly assigned over uniform intervals (i) 0 < p 0 < 0.25, (ii) 0 < p 0 < 0.5, (iii) 0.25 < p 0 < 0.75 and (iv) 0.3 < p 0 < 1.0, demonstrating the context- 
dependence of competition, (c) Osmotic competition between 30 model oleate vesicles, each with 120 nm surface diameter and each starting with a 
different degree of swelling, from maximally swelled to 5% deflated (<1> 0 randomly assigned over uniform interval 0.7 < <1) 0 < 1.0348). Any vesicle starting 
at tension state <S>o > <£[,"' = 0.853 gains lipids from its competitors. Initial fatty acid concentration inside and outside vesicles was [I] = 6.67 X 10~ 5 M. 
Environment volume for competition in all three plots was 1.04 X 10~ 14 litres. Orange crosses show agreement with equilibrated deterministic 
simulation of the model. 



appearance and eventual disappearance of competition would thus 
follow the same type of dynamics as those experimentally reported by 
Budin & Szostak 17 (Fig. 3D therein) for nervonic acid, which redis- 
tributes between vesicles. If vesicles contained a metabolism which 
synthesised phospholipid, then lasting P ;l asymmetries between vesi- 
cles could conceivably be maintained as steady states, despite a con- 
tinuous process of P exchange. However, in this work we took the 
route of not explicitly modelling phospholipid synthesis, to reduce 
the competition scenario to a materially-closed system which subse- 
quently settles to equilibrium. Under this condition, analysis is easier 
to perform. In summary, the results of this study can be interpreted as 
reflecting the competition advantage bestowed upon a vesicle by a 
membrane phospholipid fraction given that this fraction is somehow 
maintained as constant. 

The next point that deserves discussion is the the causative role of 
the direct effect in the phospholipid-driven competition simulations 
performed with our model. Our choice for function r means that a 
DOPA : OA vesicle with 10 mol% DOPA fraction will have fatty 
acids leaving the bilayer at reduced rate r(0.1) X k out = 0.9k out when 
the direct effect is maximal (d = 1). It could be argued that other 
function choices for r could reduce fatty acid off- rate even further for 
the same DOPA fraction. However, a large direct effect is not needed 
to best fit model outcomes with experimental outcomes, and would 
actually make the fit worse. Examination of Figs 4b and 4d in fact 
shows that having only the indirect effect provides the best fit to 
experimental outcomes (quantified in Supplementary Table S2). 
On the other hand, the dynamics and stoichiometry outcomes of 
Figs 4a and 4c respectively are only improved when there is a small 
reduction of k out : a small direct effect of around 0 < d < 1 for Fig. 4a 



and with d a little larger than 1 for Fig. 4c. With the maximal level of 
direct effect provided by our function r, Supplementary Table S3 
shows that the direct effect only accounts for around 20% of the total 
vesicle surface growth. Therefore, we should conclude that in our 
kinetic treatment of vesicle competition, the indirect effect is the 
main mechanism driving vesicle growth dynamics. 

One curiosity in the results (both in vitro and in silico) is how 
DOPA : OA (p 0 = 0.1) vesicles grow continually as more OA vesicles 
are added (Fig. 4c). This is unintuitive, since the growth of the 
DOPA : OA vesicles should imply a dilution of their phospholipid 
content, which would seemingly reduce the indirect and direct 
effects, thus giving a negative feedback to eventually curb the 
DOPA : OA growth profile. The reason why our model reproduces 
this continuous growth result has to do with the mathematics under- 
lying the kinetic modelling. In the limit of infinite L f , lipids in the 
membranes of our model DOPA: OA vesicles, the inside/outside 
lipid concentration required to sustain them at equilibrium (given 
by function/, as defined by equation (13) later in the Methods sec- 
tion) tends to, but crucially never actually reaches, the CVC of pure 
oleic acid: 



lim f - 

L u -»oo 



1k 0l 



k in a L 



(10) 



This is true, even if a model DOPA : OA vesicle contains just one 
single phospholipid in the membrane. Now, as more OA vesicles are 
mixed with the DOPA : OA vesicles, the population becomes increas- 
ingly dominated by OA vesicles and the lipid monomer concentra- 
tion in the environment subsequently rises toward [L]° A . As this 



SCIENTIFIC REPORTS | 4:5675 | DOI: 1 0. 1 038/srep05675 



8 





Figure 6 | Wider exploration of three different vesicle competition scenarios. Relative surface growths of two vesicle populations is explored in a broader 
context for three different competition scenarios detailed by the key. All model vesicles are considered to be 120 nm diameter. Additionally, the 
DOPA : OA vesicles are considered 5% deflated upon extrusion, and only have the indirect effect present (d = 0). (a) Phospholipid-driven competition. 
Population 1, a fixed population of vesicles with initial DOPA phospholipid fraction p^ pl =0.1, is mixed 1 : R with population 2, whose vesicles have 
initial DOPA fraction ff$ fl . (b) Osmotically-driven competition. Population 1, a fixed population of initially swelled 5>{j D '' 1 =0.85 vesicles, is mixed 1 : R 
with population 2, whose vesicles are initially swelled by O^ 2 . (c) Phospholipid versus osmotically-driven competition. Vesicles with initial DOPA 
fraction p{J° pl are mixed 1 : 1 with pure oleate vesicles swelled by O^"* 2 . Blue box outlines on the 3d plots highlight when the relative surface growth is 1. 



happens, equation (10) implies that the DOPA : OA vesicles will be 
absorbing more and more L lipids, in order to grow to a size that is in 
equilibrium with the external lipid monomer concentration. The 
growth of our model DOPA : OA vesicles is thus halted only by the 
number of lipids in the system being limited to L t and not by dilution 
of the membrane phospholipid fraction. In our kinetics model, con- 
tinuous growth happens with or without the direct effect present. 

A final point worth highlighting is that when the lipid uptake 
function u given in equation (2) is not conditional, as we assumed, 
but simply 



u(0J)= e Wi-l) (11) 



for all membrane states (which denotes that even flaccid vesicles have 
differential rates of lipid uptake), then, quite interestingly, the con- 
tinuous DOPA : OA growth effect cannot be reproduced. If we use 
definition (11) for function u in function / (13), and call the new 
function f„ a it can be shown that 
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r 2k out 
Itm f„ c = ■ 



*p(i)>M 



(12) 



meaning that the DOPA : OA vesicles do not show the same contin- 
ued growth as the lipid monomer concentration in the environment 
rises toward [L]^- Rather, the DOPA: OA have much slower 
growth, and they even have a finite stable size when the outside lipid 
monomer concentration is exactly Thus, to best reproduce 

experimental outcomes, a crucial part of our lipid uptake kinetics 
was to accelerate lipid uptake only in osmotically stressed vesicle 
states, not in flaccid ones. This is a new addition to our general kinetic 
model, introduced in this paper. 

This work is a step forward in the development of semi-realistic, 
coarse-grained descriptions of phenomena that, in reality, are extre- 
mely complex. Self-assembly processes involving heterogeneous 
component mixtures and the formation of dynamic supramolecular 
structures that could hypothetically lead to biologically relevant 
forms of material organization, like protocells 28 ' 29 , constitute a tre- 
mendous challenge, indeed, both for experimental and theoretical 
'systems chemistry' research 30 and for synthetic biology 7,31 . In par- 
ticular, the connection between basic metabolic reaction networks 
and membrane dynamics (including stationary growth and division 
cycles 32 ) needs to be explored much more extensively, since it is one 
of the key aspects to establish a plausible route from physics and 
chemistry towards biological phenomenology. 

Methods 

Fast computation of competition equilibrium. Here we provide a general numerical 
approach to solving the equilibrium configuration of a possibly heterogeneous 
population of vesicles competing for a limited supply of lipid. These vesicles may be 
osmotically swelled, laden with phospholipid, or a mixture of both, and can be 
arbitrary in number. The method allows the lipid uptake and release functions u and r 
to take arbitrary forms, subject to some requirements detailed below. 

We start by defining a function/:L /( — » [L] eq , like equation (7), which gives the 
inside/outside lipid monomer concentration [L] eq necessary to maintain a particular 
vesicle Vj at equilibrium, given that this vesicle has a specific number of lipids/ 
phospholipids in the membrane, and a specific volume: 



kin L'^+P^p u(Oj) 



(13) 



The inverse of this function yields useful information: it is the mapping of [L] eq to the 
number of lipids which must exist in the membrane of a particular vesicle, in order for 
that vesicle to be at equilibrium. 

However, due to the difficulty in isolating L fl from the potentially non-linear 
functions u and r, in most cases the inverse mapping is not possible to write in closed 
form. Nevertheless, if uptake and release functions u and r make function/(i) one-to- 
one, (ii) onto and (iii) continuous, then it follows that the inverse mapping is a 
function/ -1 , which can be numerically calculated for vesicle V, by using/ and binary 
searching for an L fl which satisfies: 



(14) 



using appropriate search bounds {normally: L"" M — 0, L™ ax — L t ). 

Crucially, having a means to calculate/ -1 gives a way of determining the total 
number of lipids existing in all equilibrated vesicle membranes, given that the inside/ 
outside lipid monomer concentration in the heterogeneous vesicle mixture is [L] eq . 
For each [L] eq , we know that each vesicle has a unique number of membrane lipids L fl , 
because/ -1 is itself one-to-one. This means that a certain [L] eq can only admit one 
single equilibrium configuration of vesicles, not multiple equilibrium configurations, 
and this lack of ambiguity is a desirable property for the method. 

The lipid monomer concentration [L]* inside/outside all vesicles in this single 
equilibrium configuration can be found by making use of the lipid conservation 
principle in equation (5): 



Er 



(15) 



[LLN A a-L(=0 



That is, at [L]*, the lipid making up the membranes of all equilibrated vesicles, plus 
the lipid monomer inside and outside the vesicles is equal to the total lipid in the 
system L t set by the initial condition. Expression (15) can also be solved by binary 



search of [L] eq between appropriate bounds, normally [L] ™" — max (f [L^ — 0) ) over 
all vesicles, [L]'" q x — min(J[L fl = L t )) over all vesicles. Finally, knowing [L] * allows to 
fully reconstruct the final sizes of all vesicles at equilibrium by substituting [L] eq = 
[I]* into equation (14) for each vesicle Vj. 

In the equilibrated population, some vesicles will have grown larger in surface area 
at the expense of others which will have shrunk. When a population of vesicles has 
competed for lipid via phospholipid -driven competition, the 'tipping point' is the 
critical number of membrane phospholipids P c *' 1 separating those vesicles which have 
lost lipid from those which have gained lipid, and is found by: 



\LV,p.,a, 



0 



(16) 



where S fL — S^, again solvable by binary searching, this time in the range 

0 <Pf t < ^p, (from a pure lipid membrane to a pure phospholipid membrane). 

Expression (16) amounts to asking how many phospholipids a hypothetical vesicle 
would require in order not to grow in surface area when the lipid monomer con- 
centration has stabilised at [L]*. Likewise, the number of phospholipids required to 
achieve any arbitrary surface area growth can be found by setting S» to the value 
desired. 

The critical phospholipid number can be stated more usefully as the critical 
phospholipid molecular fraction 



iS^ + Pf (« l -*p) 



(17) 



a vesicle has in the initial condition, a time when all vesicles have a surface of S° . For 
osmotically- driven competition, the critical volume separating shrinking vesicles 
from growing vesicles is found by searching expression (16) for vesicle volume 
instead: 



£T'< =n[f- 1 ([L]*,p„,n) - 2S " p " ap =o 



(18) 



where 5 /( — S fl . This may be alternatively stated as the critical O in the initial condition: 

S?, 



{/367t(f2" i ') 2 



(19) 



If no sign change results when evaluating the functions given in equations (14-18) at 
the upper and lower search bounds, then the respective equation cannot be solved by 
this numerical bisection approach. Otherwise typically 30 iterations of binary search 
were used to converge to an accurate answer. 
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